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Abstract. Ito equations are derived for simple stochastic cellular 
automaton with parameters and compared with results obtained 
from the histogram method. Good agreement for various param- 
eters supports wide applicability of the Ito equation as a macro- 
scopic model. 
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I. Introduction 

The Ito equation describes evolution of a stochastic diffusion Markov 
process (of order 1). Its parameters can be nonlinear functions of the 
process. Therefore, the Ito equation is a good candidate for a nonlinear 
model of phenomenon which manifests non-regular, random behaviour. 
Quite often observed variables have macroscopic character, and hence 
the Ito equation can be considered as a macroscopic model of the com- 
plex system, in which microscopic collective interactions have been av- 
eraged to an adequate form of terms in the equation. 

Various complex processes, geophysical, biological, economical etc., 
in spite of apparent random appearance, can lead to some regular be- 
haviours or patterns. Registered time series are investigated by differ- 
ent, less or more sophisticated methods, which help to find these de- 
pendencies. Some interesting regularities may constitute of new macro- 
scopic laws of complex systems. 

Time series analysis is well developed branch of science. In the linear 
case, well known procedures were elaborated (ARMA etc.). However, 
effective methods of construction of nonlinear models from time series 
data still are not very satisfactory and require further improvements. 
Using the Ito equation may provide some progress in the field. 

First attempts of determination of the Ito equation were proposed 
by Haken and Borland (Haken 1988, Borland and Haken 1992a, b, 
1993, Borland 1996)). Their procedure (SEQUIN) uses the knowledge 
of some moments of the joined distribution function, which form con- 
straints for the Maximum Information Principle. The SEQUIN method 
works sufficiently well for the case of a weak noise. However, it was 
shown (Rozmarynowska 2009) that the procedure fails in the case of 
strong multiplicative noise, when long-tail distributions appear. 

l 
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A purely numerical procedure of construction of the Ito equation 
from time series data was proposed by Siegert et. al. (1998). This 
direct procedure, based on the histogram of the joined distribution 
function, always produces an output (i.e. clouds of points) It gives some 
approximation of terms in the Ito equation. It is effective for strong 
noise too, but fitting proper functions to scattered clouds of points is 
a difficult task. Moreover, it should be underlined that this method 
may be fallacious, when the time series can not be approximated by a 
diffusion Markov process. It is possible to determine the Markov order 
of a given time series (Racca et al. 2007) but there is no method to 
verify if it is diffusion process or not. The only inspection we may 
perform is a comparison of the input time series with that generated 
by the reconstructed Ito equation. 

An adequacy of approximation of a natural phenomena by nonlinear 
models is of a great importance. To touch the topic in the context 
of the Ito equation we propose as follows. First, we replace the natu- 
ral phenomenon by the domino cellular automaton, which can be fully 
monitored. Due to a complete knowledge of the three level hierarchy 
of the model, relations between the macroscopic (the Ito equation, mo- 
ment equations) and the microscopic (rules of the automaton) descrip- 
tion are clear and understood. Then, we construct the Ito equation in 
two independent ways. The standard histogram method is compared 
with results derived analytically. Derivation of the Ito equation out of 
a cellular automaton is a crucial step in this approach. To this aim a 
special stochastic cellular automaton with avalanches was introduced 
(Bialecki and Czechowski 2010b, Czechowski and Bialecki 2010a). The 
defining rules of the domino automaton were chosen to satisfy two 
opposite requirements: complex behavior and relatively simple math- 
ematical structure. Introducing efficiency parameters in the domino 
automaton leads to a big diversity of states, which covers wide range 
of avalanche sizes. The goal of this paper is an analytical derivation 
of the Ito equations for a domino cellular automaton with efficiency 
parameters and a comparison with results of the histogram method. 

In Section[2]we shortly introduce the domino cellular automaton with 
efficiency parameters and present equations describing its behaviour in 
equilibrium state. Properties of the automaton for various efficiency pa- 
rameters are analysed. In Section Owe consider fluctuations around the 
equilibrium and propose suitable approximation formulas. In Section 
H] the Ito equation is derived. The comparison of histogram procedure 
results and analytical formulas are presented. We conclude with short 
remarks in Section [5j 
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Figure 1. Domino cellular automaton at the quasi- 
equilibrium state: left - time series for density p(i), right 
- time series for avalanche sizes w(i). Three examples for 
different cases: r] = 4 (upper diagrams), 77 = 1 (middle) 
and rj = 0.25 (lower). 



2. Domino cellular automaton with efficiency 

parameters 

We present shortly rules governing the automaton and some relations 
derived in Bialecki and Czechowski (2010a, b). Consider 1-dimensional 
grid of cells. There are only two states of a cell: empty or occupied 
(by a single ball). The evolution rules are as follows. In a time step a 
single new ball is added to the system to the randomly chosen cell. 

• If the new ball hits an empty cell it becomes occupied with prob- 
ability u, or the ball rebounds from the grid with probability 
(l-i/). 

• If the new ball hits an occupied cell it rebounds with probability 
(1 — (j,), or with probability ji it knocks out balls from the cell 
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and other adjacent occupied cells from the whole cluster (they 
become empty - an avalanche is triggered). 

An important macroscopic variable, which describes the state of the 
automaton is the density p i.e. the rate of number of occupied cells 
to the grid size. We monitor two time series generated by the cellular 
automaton: the avalanche sizes w{i) and the density p(i) computed 
after each avalanche. Figure Q] presents three examples of such time 
series for three chosen values of the ratio p/u = 4,1, and 0.25. An 
increase of the average density, of the average avalanche size and of the 
range of fluctuations of p with decreasing pjv is evident. 

A state of the automaton can be characterized by the cluster size 
distribution n i; i.e. by numbers of clusters of size i (normalized by 
the size of the grid). The following relations, which are exact in the 
equilibrium state, were derived: 

1 2 1 
(1) m = (l-p) n+-n? 

2 / n?\ 
i / / o \ 



(3)n, = ^ 2^ 1-^ +n ?V^ for i > 3. 



Here rj is the rate /i/z/, n is the number of all clusters, n = Yjci n ^ an< i 
p is the density of occupied cells, p = J2i>i ^ n i- The is the number 



,0 _ 2n 



of single empty cells and is given by the following formula — -, — 

The following simple balance equations for moments of n« were also 
derived: 

(4) 2m + (1 + rfimi = 1, 

(5) mi + 77777-2 = 1, 

where: mo = n, m\ = p and 1712 = Xli>i n ^ 2 - The average cluster 
size < i >, the average avalanche size < w > and the average square 
deviation from the average cluster size < T > (i.e. the analogue of 
temperature) are expressed by these moments: 



(6) <i> 

(7) < w > 



mi 2p 



m l-(l + 77)p 
m 2 l-p 



nil VP 

(8) <T> = y^Uiji— < i >) 2 = p(< w > - < i >). 



i>l 
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Pe 


< w > 




analytical simulation 


analytical simulation 


rj = 4 


0.1501 0.1498 


1.4147 1.4159 


rj — 1 


0.3075 0.3068 


2.2520 2.2421 


77 = 1/4 


0.4815 0.4801 


4.3074 4.3382 



Table 1. Comparison of values of equilibrium density 
p e and average avalanche size < w > calculated from 
equations with those obtained from simulations. 
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Figure 2. Dependences of equilibrium variables: p e , n, 
< i >, < w > and < T > on the ratio of efficiency 
parameters rj = pjv. 



The equilibrium value of the density p e can be found by numerical 
solution of implicit algebraic equation 

(9) p = ^2in i} 

where n\ are given by equations ([I])-®. 

Figure [2] shows dependencies of equilibrium variables: p e , n, < i >, 
< w > and < T >, on the ratio rj. As it was expected, they decrease 
with increasing t], because greater probability p (relevant to v) of 
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rj — 4 1] = 1 rj — 0.25 



a x (p) 
a 2 (p) 

a 4 (p) 



0.9528p 
0.6449p + 0.0715 
2.4095p- 0.1846 
0.5163p + 0.0989 
0.5139p + 0.0994 



0.8365p 
0.7032p + 0.1223 
1.4316p- 0.0484 
0.5651p + 0.2273 
0.5569p + 0.2328 



0.7348p 
0.5669p + 0.1470 
1.3742p - 0.0546 
0.4908p + 0.3761 
0.4772p + 0.4054 



Table 2. Numerical values of functions a k for expres- 
sions (fIUD-(fl2l. 



triggering the avalanche decrease the number of occupied cells on the 
grid. Only the total number of clusters n, has a maximum for rj around 
1/4. Tabled] display a comparison of values of p e and < w > calculated 
from equations with respective values taken from simulation. 

3. Formulas for fluctuations around the equilibrium 

For deviations from the equilibrium state p e formulas from Section 
[2] are not valid. Therefore, we use the approximation introduced in 
our parallel paper (Czechowski and Bialecki, 2010b). It was based 
on the assumption that the following geometric-like form for rii(p) is 
maintained: 

(10) m(p) = (l-p) 2 a!(p) = (l-p)Vp, 

(11) n k (p) = n fc _i(p)a fc (p) for 4 > k > 2, 

(12) n k (p) = n k ^{ P ){a 5 {p)) k ^ for k > 5, 

where a k (p), k = 1, . . . , 5 are linear functions of p. The forms of the 
formulas (|T0|) - (|T2|) extend the 1-D percolation geometric cluster size 
distribution, n k = p h (l — p) 2 . Numerical simulations of the domino 
automaton confirm a validity of this approximation. 

Linear functions a k {p), k = 2, . . . , 5, are found from the Taylor ex- 
pansion of ratios: 

(13) n fe (p)/n fc _i(p) 

around p e to the order 1. Using the procedure presented in Czechowski 
and Bialecki (2010b) numerical formulas for a k (p), k = 1,...,5, for 
chosen values of the parameter r\ were found and they are presented in 
Table [2J As a check for these expressions, we we plot the first moment 
of n k (p) (which should be equal to the density p) as a function of p for 
chosen parameters r]. Figure [3] shows satisfactory fit in the full range 
of variability of p for all three values of r\. 

4. Derivation of the Ito equation 

First of all, we should decide which time series, of the density p{i) or 
of the avalanche size w(i), may be treated as a realization of the Markov 
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Figure 3. Sum J2 i>1 in i computed by using formulas 
(IT0"|) - (IT2"|) versus p for three cases: rj — 4 (left graph), 
rj — 1 (middle), rj = 0.25 (right). 

diffusion process. Because the avalanche time series is more complex 
and not Markovian of order 1, we use the time series for density. An 
easy check ensures that automaton rules provide the required 1st order 
Markov property for p(i). We emphasize that the following convention 
is assumed: the state of the grid (and the size of the avalanche) is mon- 
itored after each avalanche only. In this way we avoid less interesting 
stairs-like increase of p(i) in periods between avalanches. 
In order to derive the Ito equation 

(14) dp = a(p)dt + y/b(pjdW(t), 

we need the transition probability function P(p + Ap, t + At \ p,t) 
(see Risken 1996). Functions a(p) and b(p) are known to be the drift 
and diffusion coefficients respectively, and W(t) is the Wiener process. 
Following our parallel paper (Czechowski and Bialecki, 2010b, in which 
the case i] = 1 was considered) the probability of effective gain EG{k) 
(an increase of p(i) by k occupied cells) and effective loss EL(k) are 
given by the formulas: 

/ ^ 

EG(k) = p( Pl + — ,i + l\ Pi ,i 

(15) = J2 Hl-p)) s ^pw s - k {p) for k>0, 

s=k+l 

( k 

EL{k) = Pl Pi -—,i + l\ Pi ,i 

(16) = £(1/(1 - p)) s - k iipw s {p) for k>l, 

s=k 

where w s = ^ is the probability that the occupied cell is a part of the 
cluster of size s (which also corresponds to the probability of avalanche 
of size s) and N is the grid size. 

Balls rebounded off the grid do not trigger avalanches, they set the 
relative efficiency for triggering an avalanche with respect the efficiency 
of occupation of an empty cell. Hence, for time series analysis, one can 
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consider only active (not rebounded) balls introducing appropriate cor- 
rective coefficient. It follows, in expressions ( II 5p and ( I16p parameters 
p and v should be replaced by effective parameters: 



(17) p e 

(18) v P 



V 1 



pp + v{l-p) rip + l-p 

As a result, probabilities EG(k) and EL(k) satisfy the normalization 
condition. 

The drift and diffusion forces in the Ito equation correspond (see 
Risken 1996) to the first and the second moment of the transition 
probability: 



(19) a(p) oc k (EG(k) — EL(k)) , 

k>l 

(20) b{p) oc J2 k2 ( EG (k) + EL{k)) . 

k>l 

Tedious transformations on sums (see Appendix in Czechowski and 
Bialecki, 2010b) lead to expressions for a(p) and b(p) in terms of the 
first, the second and the third moment of rik(p): 

(21) a {p) oc ^^-Ig**, 

(22) b(p) oc ^(l-p)(l + ve(l-p)) 

^i_ I/g (i_ p ))L* n » + p L* n »" 

FV ey F >> k>l r k>\ 

Note, the dependence on v e is not essential. The above formulas can 
be also expressed in terms of p e ; in fact they depends on r\ and p 
only. The second and the third moments itself can be calculated using 
approximate formulas ( TT0|) -( fT2|) from Section [3j 



5. Conclusions 

Introduction of efficiency parameters p and v essentially enriched 
the automaton. By changing the ratio rj = pjv the domino cellular 
automaton can attain a wide range of values for important variables: 
p e , n, < i >, < w > and < T > (see Figure |2]), which character- 
ize a macroscopic behaviour. The time evolution of the model can be 
described by the Ito equation, therefore, we derive analytically a sat- 
isfactory approximation of Ito terms a(p) and b(p) as functions of the 
ratio 7] = pjv. For three chosen values of 77 analytical and simulation 
results were compared. Figure @] shows a good fit for these three cases. 
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Figure 4. Illustration of functions a(p) and b(p) in Ito 
equations for three values of parameters: 77 = 4 (upper 
diagrams), 77 = 1 (middle), 77 = 0.25 (lower). Dots rep- 
resent values reconstructed by the histogram procedure 
from time series generated by the domino cellular au- 
tomaton (see Figured]), lines - calculated from formulas 
(12"T|) and (1221 using approximation (flu]) - (fi2l and param- 
eters presented in Table [2J 



Thus, we have a simple, fully elaborated model which can manifest a 
wide range of behaviours. 

Moreover, we show that the histogram method of reconstruction of 
Ito equation from time series works well in the case of domino cellu- 
lar automaton. The automaton plays a role of a complex system, in 
which energy is slowly accumulated and rapidly released in avalanches. 
Therefore, this is possible that the histogram procedure can offer an 
adequate method of nonlinear modelling of similar natural phenomena 
basing on the time series data. The method fills a gap between linear 
stochastic models (ARMA etc.) and nonlinear deterministic models 
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(Takens method, see Takens 1981) because it is both stochastic and 
nonlinear. 
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